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We compare Evolutionary Algorithms with Minima Hopping for global optimization in the field 
of cluster structure prediction. We introduce a new average offspring recombination operator and 
compare it with previously used operators. Minima Hopping is improved with a softening method 
and a stronger feedback mechanism. Test systems are atomic clusters with Lennard- Jones interaction 
as well as silicon and gold clusters described by force fields. The improved Minima Hopping is found 
to be well-suited to all these homoatomic problems. The evolutionary algorithm is more efficient for 
systems with compact and symmetric ground states, including LJ150, but it fails for systems with 
very complex energy landscapes and asymmetric ground states, such as LJ75 and silicon clusters 
with more than 30 atoms. Both successes and failures of the evolutionary algorithm suggest ways 
for its improvement. 



INTRODUCTION 

To find the structural ground state of a cluster is a 
non-trivial global optimization task. One has to find the 
global minimum of the potential energy surface which is 
a function of all the atomic coordinates. Even for a rel- 
atively small cluster of 30 atoms the configuration space 
has already 90 dimensions. Because knowing the struc- 
ture is a prerequisite for the study of all other physical 
and chemical properties the problem is of great impor- 
tance and many algorithms have been developed to solve 
this global optimization problem. We compare Evolu- 
tionary Algorithms which have successfully been used in 
many diverse fields with the Minima Hopping [l[ method. 

For the prediction of the ground state structure of crys- 
tals, the USPEX (Universal Structure Predictor: Evolu- 
tionary Xtallography) 0, Q method has turned out to 
be extremely powerful and has already allowed material 
scientists to find interesting and unexpected new crystal 
structures [1, 0, @, 0, Hi- Recently Evolutionary Algo- 
rithms have also been successfully used to predict sur- 
face phenomena such as steps on silicon crystals [tl,[ll|. A 
widespread application of global optimization methods is 
the prediction of the structure of various clusters, In this 
field the majority of the work has been done with genetic 
or evolutionary methods as well [II [II, [ll, [13 ■ We note 
that different evolutionary algorithms developed for var- 
ious types of structure prediction problems (molecules, 
clusters, crystals) have significant differences. Even for 
the same type of problem (e.g. crystal structure predic- 
tion) the previously proposed algorithms are very differ- 
ent in their construction and performance. 

The minima hopping method has been successfully ap- 
plied to benchmark systems [H as well as to silicon clus- 
ters and AFM tips 0. 

The presented evolutionary algorithm is similar to and 



inspired by the structure prediction algorithm USPEX. 
But since we work with non-periodic systems without a 
crystal lattice, some modifications of the original method 
were required. The version of minima hopping we are us- 
ing is based on an improvement of the two key features 
of the original minima hopping method [l[. The feed- 
back mechanism is enhanced and the Bell-Evans-Polanyi 
principle flil ] is exploited in a more efficient way by mov- 
ing preferentially along soft directions in the molecular 
dynamics (MD) part of the minima hopping algorithm. 

Our comparison of Minima Hopping and Evolution- 
ary Algorithms is based on Lennard-Jones systems, es- 
pecially the cluster with 55 atoms, which is an example 
of an easy one-funnel structure, and the 38-atom sys- 
tem, which is known to have a complicated double-funnel 
structure. We also apply the algorithms to more realis- 
tic systems, namely silicon clusters described by a force 
field and gold clusters described by an embedded atom 
potential. 

This paper is structured as follows: We first introduce 
the evolutionary method used, after a quick introduction 
to Minima Hopping and its modifications we present the 
results section containing a comparison between Minima 
Hopping and the Evolutionary Algorithm and finally we 
also test different aspects of the EA and MH. 



THE EVOLUTIONARY ALGORITHM 

Evolutionary Algorithms (EA) implement a very sim- 
ple model of biological evolution. They work on a set of 
samples — a population — which is gradually improved 
by selection and reproduction of fit members of the pop- 
ulation — individuals. Each individual is a solution can- 
didate. A single iteration step leads from a population 
to the next and is called a generation. The algorithm 
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optimizes the fitness function — in our context the neg- 
ative energy of the configuration. The operators applied 
to the population to obtain the next generation are the 
heart of the algorithm as they determine its quality and 
properties. 

In contrast to the original simple genetic algorithm, 
modern applications in the field of chemical structure 
prediction all use real value encoding instead of binary 
strings and phenotypical operators acting directly in real 
space instead of gene modification lj, ll7( • Also state of 
the art is the application of local optimization to each 
individual thus reducing the search space to basin bot- 
toms. Local optimization is done using standard tech- 
niques such as steepest descent and conjugate gradient 
methods. 

FIG. Q] presents an overview of the Evolutionary Algo- 
rithm used. It starts with a randomly initialized popu- 
lation, this is our generation 0. In each generation the 
algorithm goes through three steps: Selection, Applica- 
tion of Operators and Acceptance. The steps Selection + 
Operator produce new offspring and mutations and put 
them into an intermediate population. In Acceptance the 
next generation is selected out of all available offspring 
together with individuals from the old population. 

— ~ Population ~ 1 
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FIG. 1: Evolutionary Algorithm: A population i is 
evolved into the next generation i + 1 using mutation 
and recombination operators. E: elitism (number of 
individuals kept from the old generation). O: total 
number of offspring produced. P: population size. 

The algorithm possesses many tunable parameters. 
Most important values are populationsize, the num- 
ber of offspring produced, the number of individu- 
als kept from the last generation (elitism) and the 
mutationrate, for a complete list see TABLE HI 

As it is often the case in evolutionary algorithms in 
this field an energy slot restriction allows only one candi- 
date per energy interval energyslot in the population. 
This method of preventing multiple copies of the same 
configuration in the population may be dangerous for it 



might reject an important candidate having almost the 
same energy as an individual already known. Using force 
fields allows to calculate energy and forces with very high 
precision since the numerical noise is extremely low. It is 
thus easily possible to identify structures by their energy. 

Recently, a different structure of EA, optimized for 
parallel machines, has been presented [18[- Instead of 
a stepwise evolution this approach handles a big pool of 
individuals which is subject to continuous application of 
operators. This is closer to a biological population with- 
out sharply defined generation gaps and it solves the load 
balancing problem in parallel implementations of EA. 



Operators 

Operators are used to evolve a population to a next 
generation. We use two different kinds of operators, 
heredity operators which take two individuals as input 
and produce a child sharing properties of both parents. 
The second kind is an operator applied to a single indi- 
vidual altering its configuration [mutation). 

The selection step determines to which individuals the 
operators are applied, it is dependent on the operator. 
For a heredity operator there are two parents selected 
whereas for mutation operators only one individual is 
chosen. Selection is done using a linear ranking scheme. 
Individuals are sorted with respect to their fitness values 
and then assigned a probability depending linearly on the 
rank i. The probability of selecting the individual with 
rank i is in this case 



P[i] = Pi - (* - 1 



Pi 



(1) 



where i is the rank, starting at 1, c the parameter cutoff 
and Pi the first selection probability determined by nor- 
malization constraint. The cut-off value is the last rank 
with a selection probability above zero, all following 
ranks are assigned zero selection probability. The same 
method is also used in USPEX method 0, S|, it has 
turned out to be more efficient than Boltzmann selec- 
tion where the selection probability follows a Boltzmann 
distribution depending on relative fitness values. The se- 
lection of the same individual serving as both parents is 
prevented. Mutation operators are applied randomly to 
the whole population. 

The first heredity operator used is the cut and 
splice (cutting-plane) method introduced by Deaven and 
Ho [n], E3. Both clusters are centered at their center 
of mass (COM). A randomly oriented plane cuts the two 
clusters apart. The new offspring cluster consists of one 
half of the first and the other half of the second cluster. 
Though able to obtain two offspring by this method we 
only produce one. The plane usually contains COM and 
its cut preserves the total number of atoms in the child 
cluster. 
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Function 


Name 


Standard Value 


Population size 


populationsize 


30 


Number of offspring produced 


offspring 


populationsize 


Number of individuals taken from former population 


elitism 


1/3 populationsize 


Last rank with selection probability > 


cutoff 


2/3 populationsize 


Relative rate of offspring produced with average method 


avgoff spring 


0.50 


Only one individual allowed within this energy interval 


energyslot 


KT 4 


Total rate of mutation 


mutationrate 


0.05 


Random walk mutations (relative to total mutations) 


mutrwalk 


0.60 


Strain mutations (relative to total mutations) 


mutstrain 


0.30 


Probability of random rotation before recombination 


raterndrot 


0.90 


Convergence criterion for force norm in local optimizer 


f nrmtol 


10 -4 



TABLE I: Parameters of the Evolutionary Algorithm 



A new way of producing offspring is implemented in the 
average offspring method. Both clusters are centered at 
COM and for each atom of the first cluster the closest 
lying atom of the other cluster is identified. The atom of 
the child is now placed randomly on the connecting line of 
the two parent atoms. The randomness of this operator is 
necessary to prevent producing a lot of identical offspring. 

Before application of either heredity operator a random 
rotation can be performed on one of the parents. The 
frequency at which this rotation is used can be adjusted 
via raterndrot. 

Mutations are introduced to keep the diversity of 
the population high and prevent premature convergence. 
Three different methods are used. The easiest of those is 
the random walk mutation where atoms are randomly 
displaced. The displacement is approximately normal 
distributed with a mean displacement in the order of the 
two-body potential equilibrium distance (bond length). 
A strain mutation applies a geometrical deformation to 
the whole cluster, inspired by mechanical stress. The de- 
formations include (anisotropic) compression and shear. 
Such strain transformations gave increased efficiency in 
the USPEX method. If the compression is high this 
method relates to the Big- Bang- Algorithm where con- 
figurations are relaxed from very high compressions [p| . 
We only apply a moderate compression. In a third type 
of mutation the cluster is cut into two pieces similar to 
the plane-cutting method. One of the pieces is rotated 
around an axis perpendicular to the cutting plane by a 
random angle. This method was introduced as twinning 
mutation [20(. 

Since our operators are able to produce configurations 
with atoms lying very close to each other we use a pre- 
relaxation method which is essentially a steepest descent 
with a very small step size. After a few steps the real 
self-adjusting steepest descent is started until a nearly 
quadratic region around the local minimum is reached. 
The final minimization is then carried out by a conjugate 
gradient method. 



If an intermediate set of offspring has been created 
the Accepting step is triggered. In this step it is decided 
whether an offspring is accepted into the new population 
or is discarded. The algorithm first accepts the best indi- 
viduals from the former population {elitism). The num- 
ber of individuals chosen that way is given by elitism. 
In a second step the individual with the worst fitness 
value is replaced by the best offspring if energy slot con- 
straints are fulfilled. This is repeated until all offspring 
are processed or the population is complete. 

The algorithm is left running until a given limit of gen- 
erations has been reached or the (known [l3| ) global min- 
imum has been found. 



MINIMA HOPPING 

Minima Hopping is a recently developed global opti- 
mization algorithm [l[ which makes use of a Bell-Evans- 
Polanyi (BEP) principle for Molecular Dynamics (MD) 
trajectories [16|. The BEP principle states that low en- 
ergy MD trajectories are more likely to enter the basin 
of a lower lying adjacent minimum than high energy tra- 
jectories. The algorithm also incorporates a history to 
repel it from previously visited regions. Using local op- 
timization and molecular dynamics simulation it jumps 
between basins of attraction. The kinetic energy is kept 
as low as possible to escape the local minimum but is in- 
creased if this minimum has already been visited before. 

Minima Hopping works with two self-adapting param- 
eters, the kinetic energy of a MD escape step ekin and 
an acceptance threshold ediff for new minima to in- 
troduce a further downward preference. Starting from a 
local minimum a MD escape trial is started with kinetic 
energy ekin. After a few steps it is stopped and the con- 
figuration locally optimized again. If escaped the new 
minimum is only accepted when the new energy lies at 
maximum ediff higher than that of the previous mini- 
mum. 
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Minima Hopping is adjusted by tuning five feedback 
parameters: a\ decreases edif f when a new minimum 
is accepted whereas a 2 increases the threshold on rejec- 
tion, Pi increases the kinetic energy when a MD escape 
trial fails, /3 2 increases ekin when the new minimum is 
already known and (3% decreases the kinetic energy if the 
minimum is unknown. Each visited minimum is added to 
a history list and marked as known. The algorithm cur- 
rently uses the energy value to identify different minima. 
For more details we refer to the original paper [l[ . 

Minima Hopping was used with the standard parame- 
ter set presented in the original paper: ai = 1/1.05, a 2 = 
1.05 and /3i,/? 2 = 1.05, #5 = 1/1.05. 

The algorithm is efficient since it inhibits revisiting 
the same configuration many times which is likely to be 
the case in thermodynamically inspired methods such as 
Simulated Annealing. 

We present two modifications of the original algorithm: 
The initial velocity vector of a MD escape trial is moved 
towards a direction with low curvature (softening) and a 
stronger feedback mechanism is used. 

Softening. MD escape trials in the MH algorithm 
need an initial velocity distribution which is then rescaled 
to fit the desired kinetic energy. The velocities are ran- 
domly directed for each atom with Gaussian distributed 
magnitudes. Regardless of the actual distribution chosen 
it has proved very useful to use softening, to choose ve- 
locities along low-curvature directions. In this way one 
can typically find MD trajectories with a relatively small 
energy that cross rapidly into another basin of attraction. 
In the original MH method low kinetic energy trajectories 
could only be obtained by using large values for mdmin 
which results in long trajectories. A direction of low cur- 
vature is found using a modified iterative dimer method 
which only us es g radients, no second derivatives need to 
be calculated [21| . 

Starting at a local minimum x with an escape direction 
N the method calculates a second point y = x + rfN 
at a distance d along the escape direction. The forces 
are evaluated at y and the point is moved along a force 
component F 1 - perpendicular to N: 

p-L = F - (F ■ N)N 



y = y 



aF J 



N' = 



y x 



After a few steps the iteration is stopped before a locally 
optimal lowest curvature mode is found. Initial velocities 
for the MD escape are then chosen along the final escape 
direction N. 

If the softening procedure is executed until it converges 
the performance drops again. It is important not to 



overdo softening. Always escaping into the same soft 
mode direction of a given minimum reduces the possibil- 
ities of different escape directions and therefore weakens 
the method. A good indicator was the mean kinetic en- 
ergy during a run. For a few softening iterations the 
value decreases whereas it starts to increase again at a 
certain number of softening iterations. We set the iter- 
ation count to the value where the mean kinetic energy 
was minimal. 

Typically 40 iterations are done with a step size a and 
a dimer length of d = 0.01. 

Enhanced Feedback. In original MH ekin is increased 
by a factor /3 2 if the current minimum has already been 
visited before — regardless of the number of previous 
visits. An enhanced feedback method uses a value of /3 2 
depending on the previous visits according to 



ft =/3 2 ° x (1 + clogAO 



(2) 



where /3 2 i s the original value of 1.05 and N the number 
of previous visits to this minimum. The parameter c has 
been set to 0.1 after tests on bigger Lennard- Jones clus- 
ters and gold systems. This feedback mechanism reacts 
slightly stronger if the minimum is visited many times. 
If the system has only one energy funnel this enhanced 
feedback can even be slightly disadvantageous since it 
increases the kinetic energy too much and thus weak- 
ens the BEP effect of MD. The increased feedback mech- 
anism improves the efficiency however considerably for 
large systems where the system can be trapped in huge 
structural funnels. If a cluster has for instance both low 
energy icosahedral and decahedral structures it takes a 
very long time for the MH algorithm without enhanced 
feedback to switch from one structure to the other. 

Parallel Minima Hopping. Parallelizing Minima Hop- 
ping is straightforward. On each processor an own MH 
process is started, all sharing the same history list. Only 
energy values of visited minima have to be shared. Due 
to the feedback mechanism and the common history list 
overlap of search areas is penalized in this parallel setup 
and running on several processors can thus yield an al- 
most linear speedup in runtime. A further effect can 
be exploited in parallel runs. A single run might eas- 
ily get trapped in a metastable funnel with a high escape 
time but the probability of all processors getting trapped 
in this local minimum is exponentially reduced with the 
number of processors running in parallel. The total ex- 
pected runtime until success is therefore less influenced 
by the long escape times of relatively stable local minima 
in parallel runs. On the other hand there is a minimal 
number of local minima that have to be visited before the 
global optimum can be found. This minimal number of 
hops renders the use of too many processors less efficient 
again thus leading to an optimal number of processors 
depending on the structure of the PES. The idea of par- 
allel sampling is known to have a positive effect 22 1. 
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COMPUTATIONAL EXPERIMENTS ON 
CLUSTERS 

Atomic clusters are an ensemble of bound atoms, big- 
ger than a dimer but smaller than bulk matter. They 
show interesting properties in the transition region of 
single atoms to bulk matter. From a global optimization 
point of view they are complicated multi-dimensional sys- 
tems usually difficult to optimize since they contain a lot 
of local minima. Therefore they arc of interest to test 
the capabilities of an optimization algorithm. 

A simple model of chemical interaction of two non- 
charged atoms is the Lennard- Jones potential. The po- 
tential well depth and the equilibrium distance are the 
only parameters, both are set to 1. 

Such models represent rare-gas clusters reasonably 
well. Lennard- Jones clusters are thoroughly studied 
model systems with well-known global minima up to 1000 
atoms. Some of those clusters have a non-trivial multi- 
funnel potential energy surface (PES) where the global 
minimum is not easily found. The use of the Lennard- 
Jones two-body potential is fast compared to more accu- 
rate potentials or even Density Functional Theory (DFT) 
calculations. For these reasons they are well-suited to 
test global optimization methods. Our algorithms are 
applied to Lennard- Jones systems consisting of 38, 55, 
75, 100 and 150 atoms. All clusters chosen show a differ- 
ent aspect, LJ 55 is a very easy system, LJ 38 has a double 
funnel structure, LJ75 is a very hard double-funnel sys- 
tem and LJ100 and LJ150 are examples of bigger clusters. 
Additionally, we perform single runs on LJ39, LJ74 and 
LJg 8 to find the ground state motif dependence of the 
algorithms. For these additional clusters we have pairs 
of similar-sized clusters with very different ground state 
geometries. 

Silicon clusters are of more practical interest as silicon 
is so widely used in research and technology. To obtain 
realistic results usable in those fields it would be neces- 
sary to calculate energies using at least DFT methods. 
But those methods are consuming a lot of CPU time 
and one should highly optimize the algorithms applying 
DFT. Especially in the field of global optimization it is 
of importance to use an algorithm which is efficient to 
save computing time. To investigate the behavior of the 
presented algorithm we use only a force field to evaluate 
the energy of a configuration. Silicon systems are chosen 
as model systems since they possess directed bonds and 
show frustrated behavior. They have non-trivial min- 
imum structures which are in general neither compact 
nor spherically symmetric. Silicon systems containing 18, 
22, 30 and 60 atoms are explored. The force field used is 
Bazant's Environment- Dependent Interatomic Potential 
(EDIP) [H, [H, [H] . This force field has been chosen 
because it tends to elongated minimal structures which 
are not spherical (FIG. [2]). There exist other force fields 
preferring spherical ground states. But the performance 




FIG. 2: Si 30 ground state with EDIP. 



with spherical ground states is already investigated in 
Lennard-Jones and gold systems. Such highly elongated 
structures might not be very realistic as more accurate 
DFT calculation suggest only slightly elongated ground 
state configurations. 

In gold systems different configurations can have very 
similar energy and the global optimum is often only very 
slightly lower than the next-best solution. This leads to 
a situation where the global geometry of the gold cluster 
can change completely by adding only a single atom. The 
energies of gold clusters are calculated using an empirical 
many-body potential by Rosato, Guillop and Legrand 
(RGL) f27| . Gold clusters with 28, 76 and 102 atoms are 
investigated. Additionally, AU79 and Auioi are added 
to the test set in order to have different ground state 
geometries at similar sizes. 

To compare different algorithms we measure two quan- 
tities, the total number of calls to the function calculating 
energy and forces and the total number of local geometry 
optimizations. The local optimization is the most expen- 
sive step in both algorithms and the speed is mainly de- 
termined by its number. To obtain meaningful numbers 
we perform between 20 and 100 runs on each problem. 

The numerical convergence criterion for local optimiza- 
tions is ||F|| < 10~ 4 for the total euclidean force norm 
which leads to an energy precision of almost ~ 10~ 8 . 



RESULTS 

We gathered results for total performance on the inves- 
tigated systems. Additionally, we also did comparative 
runs concerning the new modifications mentioned above. 
The numerical results of the different performance runs 
can be found in TABLE HTT1 whereas special test cases are 
addressed later. The table shows the results obtained us- 
ing the best parameter set known to us. Those values are 
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usually a combination of both heredity methods and all 
mutations. Minima Hopping found the global minimum 
in all problems investigated whereas the evolutionary al- 
gorithm could not find the global minima for LJ75, LJ 98 , 
Si 30 , Si 60 , Au 79 and Aui 02 - 

General Performance 

The comparison shows a good applicability of Minima 
Hopping in all problems. The Evolutionary Algorithm 
is capable of finding most of the ground states but fails 
in some cases. Where it is successful it can even outper- 
form MH. Problematic for the EA are the non-icosahcdral 
ground states, such as LJ75, LJgs and the elongated sili- 
con clusters. 

The failure of the Evolutionary Algorithm to find the 
elongated ground state structure of silicon systems shows 
a clear advantage of Minima Hopping, it is not geome- 
try dependent as the current Evolutionary Algorithms 
are. Moving via MD it is applicable to any system for 
which forces arc available. In the past there have been ap- 
proaches inspired by genetic algorithms which used stan- 
dard crossover and in the beginning even binary strings 
as coordinate representation. Those ideas seem to be a 
bit less geometry dependent than the nowadays applied 
heredity methods. But the original simple genetic meth- 
ods have been shown to be less effective than more elabo- 
rate geometrical operators in Lennard- Jones systems [l7| . 

When comparing two clusters of similar size but with 
different ground state structures it is obvious that the 
evolutionary algorithm can in general not reproduce non- 
icosahedral ground states with the operators presented 
here. The results of this tests can be found in TABLE HT1 
This table contains results of comparative runs without 
enough statistics to compare performance. The prob- 
lem has already been identified and solved by niches to 
prevent domination of the whole population by only one 
geometry type [28}. Minima Hopping is able to find the 
non-icosahedral motifs without further modification in 
the standard configuration though with a decreased per- 
formance comparing to the icosahedral configuration of 
similar size. 

The results show that a further adaption of EA to some 
specific features of clusters is necessary to enable the EA 
to work as efficiently for clusters as they do for periodic 
solids. 

A further difference can be noted when comparing 
function calls. Minima Hopping tends to need in total 
more calls to the energy and forces function than the 
evolutionary algorithm in the cases where the EA suc- 
ceeded, at least in smaller systems. This is clearly due 
to the use of MD and softening in Minima Hopping. But 
we remark that it is not necessary to perform MD escape 
with the full accuracy. A DFT calculation can be done 
using a reduced basis set and therefore in significantly 



System 


Structure 


Success 


LJ39 


Icosahedron 


yes 


LJ38 


fee 


yes 


LJ74 


Icosahedron 


yes 


LJ75 


Marks decahedron 


no 


LJ100 


Icosahedron 


yes 


LJgs 


Tetrahedron 


no 


Au 7 6 


Icosahedron 


yes 


Au 7 g 


fee 


no 


Auioi 


Icosahedron 


yes 


Au W 2 


fee 


no 



TABLE II: Performance of the EA depends on the 
motif of the ground state configuration. It is not 
possible to find non-icosahedral ground states in 
systems bigger than 38 atoms. Minima Hopping was 
able to find all of the listed configurations. A no means 
failure to find the ground state structure within 100 000 
local optimizations. 

less computer time. On the other hand Minima Hop- 
ping never produces energetically awful candidates since 
it moves via MD. Since heredity and mutation operators 
do not necessarily generate chemically reasonable config- 
urations the local geometry optimization requires more 
force evaluations in the EA. Since in DFT applications 
the early stages of local optimization do not need high 
precision they could even be done using force fields thus 
reducing the computational demand. 

Heredity Methods 

While comparing both heredity methods used we ob- 
served that average offspring method performed better in 
systems with a compact optimal structure. The plane-cut 
method could only produce offspring as good as the other 
method by choosing a very low cutoff. When Dcaven et 
al. [l| applied this method they produced every single 
combination of two parents from a very small population 
so having actually a very low cut-off level too. With the 
lower cut-off level the plane- cut method performed well 
overall. 

The samples obtained using average offspring method 
are energetically closer to the fittest member of the pop- 
ulation where the cutting plane method samples broader 
with more diversity (FIG. [Ha/b). It is a known fact that 
decreasing diversity in the population can lead to pre- 
mature convergence. On the other hand, if samples in 
the complete energy range are allowed the method re- 
sembles too much a random search and loses efficiency. 
To observe a well evolving population it is necessary to 
have a balanced distribution of selected individuals. It 
seems therefore necessary to lower the cut-off parameter 
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Cluster 




Minima Hopping 




Evolutionary Algorithm 




Energy 


#G0 1 


calls (GO) 2 


calls (MD) 3 


runs 4 


#GO 


calls 


configuration 


runs 


L J26 


-108.31562 


96 


50 610 


8 300 


100 


56 


34 200 


10-4-6 


100 


LJ3S 


-173.92842 


1 190 


688 500 


85 930 


100 


1265 


732 900 


25-0-20 


100 


LJ55 


-279.24847 


190 


74 840 


14 700 


100 


100 


54 900 


10-3-3 


100 


LJ75 


-397.49233 


27 375 


12.3 x 10 6 


2.11 x 10 6 


20 










LJ100 


-557.03982 


5 960 


1.87 x 10 6 


417000 


42/50 


5 908 


3.89 x 10 6 


60 - 20 - 40 


35/40 


LJ150 


-893.31026 


9 490 


3.72 x 10 6 


758 000 


45/50 


7 980 


4.36 x 10 6 


60 - 20 - 40 


17/20 


Sin 


-44.91223 eV 


29 


11790 


4190 


60 


61 


31300 


10-3-4 


50 


Siig 


-74.88419 eV 


110 


23 100 


10 700 


93 


195 


110 100 


10-4-6 


40 


Si 22 


-92.48090 eV 


370 


187400 


52 100 


21 


3 300 


1.74 x 10 6 


10-4-6 


1/5 


Si30 


-126.0952 eV 


5 050 


3.98 x 10 6 


750 000 


100 










Sieo 


-253.05509 eV 


23 300 


14.1 x 10 6 


2.80 x 10 6 


44/49 










AU28 


-99.95115 eV 


54 


26 210 


3 510 


100 


87 


68 910 


10-3-6 


50 


Auto 


-279.34791 eV 


1124 


526 500 


70 800 


98 


2 680 


1.60 x 10 6 


25-5-13 


19/20 



TABLE III: Best Results. Highlighted in bold face are the problems where the EA performed better. All values are 

averaged over multiple runs, indicated in the column runs. Standard deviations are as big as the mean values. 1 
Average number of local optimizations over all runs. 2 Calls to energy function during local optimization including 
softening. 3 Calls to energy function from MD escape steps. 4 If two numbers are present some runs had to be 
stopped due to too long runtime, the first number indicates the successful runs. 5 Configuration of EA: 

populationsize — elitism — cutoff. 



Cluster 


Average Offspring 


Plane-Cut 




#GO 


configuration 


#GO 


configuration 


LJ55 


119 


10-2-8 


100 


10-3-3 


LJ38 


1265 


25-0-20 


1595 


25-0-10 


Siis 


322 


10-3-6 


195 


10-4-6 


AU28 


87 


10-3-6 


88 


10-2-6 



TABLE IV: Heredity methods in direct comparison 



in plane-cut runs to select more of the fitter individuals 
as parents. A mix of both heredity methods delivered 
the best results (FIG. He). 

The plane-cut method is better suited to a general ap- 
plication of the algorithm, it can partially solve geome- 
tries less compact than the average offspring method. 

In big clusters it proved useful to combine all operators 
available. The results, especially LJ100, were best with 
a combination of all presented methods. In general the 
mixture was at a 1 : 1 ratio or even more preferring 
plane-cuts in systems with known tendency towards non- 
spherical ground state (e.g. Sirs in TABLE ITV|) . 

We have also tested different plane-cut setups with a 
slightly modified method where a minimal distance be- 
tween the two cluster halves is enforced. This method 
performed poorly and was always weaker than all differ- 
ent methods tested. Another modification where COM 
is not enforced to lie in the plane was also considered and 
dropped since there were no improvements. 

The random rotation before recombination is necessary 



Rotation 1 


AO 


PC 


0% 


> 1000 


280 


10% 


343 


142 


50% 


174 


142 


75% 


123 


127 


90% 


121 


115 


100% 


130 


147 



TABLE V: Number of local optimizations for different 

random rotation rates before heredity operator 
application for LJ55 systems. For each setup 20 runs 
were performed. Random rotation is crucial for average 
offspring method (AO) but not for plane-cut method 
(PC). 1 Frequency of random rotation. 



for average offspring method but only of advantage for 
the plane- cut method (see TABLE [Vj) . If the system 
prefers non-spherical configurations raterndrot should 
be small. 



Parameter Tuning 

In systems with a double-funnel structure where the 
global maximum is located in the narrower funnel it 
might turn out advantageous to disable elitism, or (bet- 
ter) include in elitism only sufficiently different struc- 
tures. In our tests LJ38 performed best when filling the 
population with offspring only. But we should remark 
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(c) Combined 1:1 

FIG. 3: Candidate samples in LJ75 after 10 000 minima, parents: all energy values of candidates selected as parents, 

offspring: all produced samples, init: initial population. The global optimum is located at —397.5. (a) Average 
Offspring method with a relative cutoff of 1, (b) Plane-Cut method with a relative cutoff 1/3, (c) Combined 1:1 

with relative cutoff 2/3. 



that average offspring method has a rather preserving 
character often reproducing candidates already known. 

A drawback of the evolutionary algorithm is the need 
to tune many parameters. A solution working with good 
performance on many different systems without adjust- 
ment would definitely be of interest. Minima Hopping 
performed better in this aspect, it never needed addi- 
tional tuning, all runs were done using the same set 
of standard parameters. Using a standard set (TA- 
BLE HI for all problems resulted in performance loss of 
the EA. The problems still converged but with perfor- 
mances down to half of the tuned versions shown in TA- 
BLE Unl A possible way to overcome this limitation is to 
fix the relative rates of elitism and mutation etc. and only 
adjust populationsize to the specific problem. Another 
possibility would be an automatically self-adapting ver- 
sion which tunes the parameters during runtime. In this 
case a stable and efficient scheme of parameter adapta- 
tion would be needed which is clearly not a trivial task. 

We note that for crystals with up 30-60 atoms in the 
unit cell, the USPEX algorithm 0, [1| proved to perform 
very well with essentially a universal set of parameters 



without any parameter tuning. Evolutionary optimiza- 
tion of clusters, which are more complex systems, is more 
sensitive to parameter values. 

Modification of Minima Hopping 

Minima Hopping has been considerably improved us- 
ing softening, in all studied cases. 

The use of enhanced feedback method is advantageous 
in large or multi-funneled systems but can even have a 
negative effect in easy systems as LJ55 (TABLE IVI[) . The 
parameter c in should not be chosen too large. We 
used softening and enhanced feedback with c = 0.1 in the 
comparison runs. 

CONCLUSION 

We have tested an Evolutionary Algorithm capable of 
finding ground state structures of atomic clusters. In 
spite of the success of EA for periodic systems and on 
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System 


No Softening 


Softening 




c = 1 


c = 


c = 0.1 


c = 0.2 


LJ38 


2 062 


1217 


1190 


990 


LJ55 


320 


140 


190 


198 


LJ100 


9100 


7300 


4 700 


5 800 


LJ150 


n/a 


14 900 


9111 


11630 


Au 2 8 


167 


44 


44 


56 


Au 76 


n/a 


979 


890 


1024 



TABLE VI: Geometry optimizations with and without 
softening in Minima Hopping with different feedback 
parameters, n/a: not tested. 1 Parameter c is defined in 
eq. © 



surfaces the current Evolutionary Algorithm is overall 
less efficient than Minima Hopping in the current imple- 
mentation. It is not yet able to find global minima with 
geometrically difficult structures such as elongated sili- 
con clusters and non-icosahedral ground states without 
the concept of niches. In contrast, Minima Hopping was 
able to find all ground states. Where the EA succeeds its 
performance is comparable to or even better than that of 
MH. 

Further improvements of the EA could make it superior 
to MH for cluster optimization if the specifics of cluster 
systems are taken into account, as it was already done 
for periodic systems in the USPEX algorithm. The MH 
algorithm, on the other hand, shows a high stability and 
little need for further adaptation, at least for homoatomic 
systems. 

Minima Hopping was improved by doing escape steps 
in directions with relatively low curvature of the PES and 
by using an enhanced feedback method. 
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